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Abstract 

The increasingly used real time quantitative reverse transcription-PCR (qRT-PCR) method for gene expression analysis 
requires one or several reference gene(s) acting as normalization factor(s). In order to facilitate gene expression studies in 
sugarcane {Saccharum officinarum), a non-model plant with limited genome information, the stability of 13 candidate 
reference genes was evaluated. The geNorm, NormRnder and deltaCt methods were used for selecting stably expressed 
internal controls across different tissues and under various experimental treatments. These results revealed that, among 
these 13 candidate reference genes, GAPDH, eEF-la and elF-4n were the most stable and suitable for use as normalization 
factors across all various experimental samples. In addition, APRT could be a candidate for examining the relationship 
between gene copy number and transcript levels in sugarcane tissue samples. According to the results evaluated by 
geNorm, combining CUL and eEF-la. in hormone treatment experiments; CAC and CUL in abiotic stress tests; GAPDH, eEF-la 
and CUL in all treatment samples plus CAC, CUL, APRT and TIPS-41 in cultivar tissues as groups for normalization would lead 
to more accurate and reliable expression quantification in sugarcane. This is the first systematic validation of reference 
genes for quantification of transcript expression profiles in sugarcane. This study should provide useful information for 
selecting reference genes for more accurate quantification of gene expression in sugarcane and other plant species. 
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Introduction 

Real time quantitative reverse transcription-PCR (qRT-PCR) is 
increasingly used in gene expression analysis owing to its simple, 
reproducible and high-throughput features. qRT-PCR provides a 
useful and rapid means of understanding gene expression in living 
organisms by measuring the expression of target genes across 
different samples. In addition, qRT-PCR is a low cost and widely 
accepted method in the tracking of gene expression levels in 
genetically modified organisms (GMO) as well as in molecular 
breeding and gene mining. When performing qRT-PCR analysis, 
several factors such as sample amount, RNA integrity, cDNA 
quality, as well as the tissues or cell activities, can affect the 
quantitative measurement of gene expression [1]. Thus, in order 
to obtain a reliable analysis of gene expression by qRT-PCR, one 
or several reference genes should serve as the internal control to 
normalize and monitor the expression variation between samples 
and reactions [2-4]. The expression of these reference genes 
should remain stable under various experimental treatments and/ 
or at different stages of development and growth periods [2-4]. 
Specifically, a suitable reference gene for performing qRT-PCR 
analysis should: (i) have stable expression across all or most of the 
samples analyzed; (ii) have no association with any pseudogene, to 
avoid the amplification of non-functional gene family members; 



(iii) reflect variations in RNA quality and quantity, as well as 
cDNA synthesis and amplification; (iv) possess the stability of 
transcription that is suitable for the target gene; and (v) exhibit 
moderate expression levels (i.e. a threshold cycle Ct of 15 to 30) 
[5]. 

Housekeeping genes (HKGs) related to basal cell activities and 
cellular structure components have been historically used as 
internal controls in medical science [6] and later in plant science 
[7]. The most commonly used HKGs are 25S rRNA, GAPDH 
(glyceraldehyde-3-phosphate dehydrogenase), ACT (fi or y actin) 
and TUB (a or [S tubulin). However, using a genome wide 
approach, hundreds of genes in Arabidopsis thaliana were shown to 
outperform traditional reference genes in terms of expression 
stability throughout development and under a range of environ- 
mental conditions [8]. Among them, genes encoding a protein 
phosphatase 2A subunit, a coatomer subunit and an ubiquitin- 
conjugating enzyme were identified as novel reference genes [8] . 

Recently, a range of new reference genes have been validated 
across a set of tissues and differendy treated samples by systematic 
statistical algorithms termed geNorm [2], BestKeeper [9], 
NormFinder [1], the deltaCt method [10] and the RefFinder 
WEB-based software (http://www.leonxie.com/reference gen- 
e.php) [11]. These genes include UBQ5 (Ubiquitin5), eEF-la 
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(Eukaryotic elongation factor 1 -alpha) and eIF-4a (Eukaryotic 
initiation factor 4-alpha) in Oryza sativa; CUL (Cullin), FPGS 
(Folylpolyglutamate synthase), LUG (Leunig), MEP (Membrane 
protein) and UBCP (Ubiquitin carrier protein) in J?ea mays; CAC 
(Clathrin adaptor complex) and TIPS- 41 (Tonoplastic intrinsic 
protein41) in Brassica juncea as well as APRT (adenine phosphor- 
ibosyl transferase) in Solatium tuberosum and Setaria italica L. [12-17]. 
Zhang et al. investigated ten candidate genes in five different 
monocot plants (Brachypodium beauv, Hordeum vulgare, Sorghum bicolor, 
Triticum aestivum and mays) under infection with different viruses. 
They found that GAPDH performed well in B. beauv and that EFlcn. 
(designated eEF-la in the present study) performed well in T. 
aestivum [18]. Additionally, a number of reference genes have been 
validated in plants such as Solarium tuberosum [19], Glycine max 
[12,16,20], Solarium lycopersicum L. [21,22], Vitis vinifera [23,24], T. 
aestivum [13,25], Brassica napus [26], cereals (7^ aestivum, H. vulgare 
and Avena sativa L.) [27], Cucumis sativus Linn. [5], Nicotiana tabacum 
[28] and Phyllostachys edulis [29]. Lastly, Caus et al. identified 
endogenous reference genes in the base genome of sugarcane 
(Saccharum qfficinarum) [30], and two of them, APRT and PRR 
(Pseudo response regulator), are used in the present study to test 
the stabilities of "low copy number genes" in transcripts. 

Sugarcane (S. qfficinarumxS. spontaneum) is a widely grown sugar 
crop in the tropics and subtropics with increasing demand due to 
biofuel production and challenges with biomass production [31]. 
Identifying genes for sucrose accumulation and stress resistance, 
which requires a set of reference genes for gene expression 
normalization, could serve to increase sugarcane yield and sucrose 
content using both genetically modified strategies and molecular 
marker-assisted breeding. Sugarcane has limited available geno- 
mic information so only a few of genes, mainly GAPDH and 25S 
rRMA, have been verified as practicable reference genes [32,33]. 
However, there is expression variation of GAPDH and 25S rRMA 
under various experimental conditions [2,12,34]. 

In the present study, the 13 genes 25S rRMA, GAPDH, ACT (R- 
actin), TUB (B-tubulin), APRT, PRR, 18S rRNA, eEF-la, eIF-4a, 
CAC, TIPS-41, CUL and UBQwere selected as candidate reference 
genes for evaluation in sugarcane. This study utilized 57 sugarcane 
samples and aimed to reveal which reference genes should be used 
in experimental samples with different treatments or different 
tissues. A combination of reference genes was also introduced to 
evaluate their potential for more accurate and reliable qRT-PCR 
analysis of gene expression in sugarcane. 

Materials and Methods 

Ethics Statement 

No specific permissions were required for these locations/ 
activities. The field studies did not involve endangered or 
protected species and the specific location of this study is longitude: 
119.23E, latitude: 26.08N. 

Plant Materials, Growth Conditions and Treatments 

Five sugarcane cultivars, "ROC"20, "ROC'22, FN40, Liu- 
cheng03-182 and YC05-179, were used for tissue sample 
collection. Five tissue samples were harvested (leaf, leaf sheath, 
stem epidermis, stem pith and bud) from 7- to 8-month-old 
sugarcane grown in the field. For each cultivar, 6 plants from the 
same experimental plot were collected to provide 6 replicates. The 
leaf and leaf sheath samples were taken from the last fully 
expanded leaf (+1 leaf), while the stem samples (epidermis and 
pith) and the buds were harvested from the 6 th or 7 th stem 
internodes. All materials were cut into small pieces, wrapped in 



tinfoil, immediately snap-frozen in liquid nitrogen and then kept at 
-80°C until RNA isolation. 

Four sugarcane cultivars, "ROC'20, FN40, Liucheng03-182 
and YC05-179, were used for the stress experiments. The single 
bud node shoots that used for in vitro disease-free plantlets 
regeneration, were incubated in 50°C water with the fungicide 
carbendazim (100 mg-L Tianjin, China) for 40 min. The 
shoots were then planted in autoclaved soil before harvest for 
meristem excision, callus induction, shoot differentiation and 
rooting [35] . The plantlets were then placed in distilled water and 
kept for ten days in a tissue culture room under a constant 
temperature of 25±1°C. Different sets of plants were then 
transferred into test tubes containing water solutions along with 
different treatments, including abscisic acid (100 uM, ABA), 
methyl jasmonate (25 U.M, MeJA) and salicylic acid (5 mM, SA) 
for 6 h; hydrogen peroxide (500 mM, H 2 0 2 ), sodium chloride 
(250 mM, NaCl) and polyethylene glycol (25% w/v, PEG) for 
12 h or copper chloride (100 |jM, CuCl 2 ) and cadmium chloride 
(500 u.M, CdCl 2 ) for 24 h. Each set of samples comprised three 
seedlings as biological replicates for expression analyses. The 
plants without treatment (kept in distilled water) were harvested as 
control. All materials were wrapped up in tinfoil, immediately 
snap-frozen in liquid nitrogen, and then kept at — 80°C until RNA 
isolation. 

A total of 57 samples (56 plus the control), including 25 tissue 
samples and 31 treated samples (the H 2 0 2 treated sample of 
YC05-179 was absent) exposed to various stress treatments 
involving four different cultivars, were employed in the experi- 
ments for evaluation of candidate reference genes in sugarcane. 

RNA Isolation, DNase Treatment and cDNA Synthesis 

TRIZOL reagent (Invitrogen Co., Carlsbad, CA, USA) was 
used in the RNA isolation of 25 tissue samples from cultivars 
"ROC'20, FN40, Liucheng03-182 and YC05-179 following the 
manufacturer's instructions. These RNA samples were treated 
with RNase-free DNasel (Promega, Fitchburg, WI, USA) before 
used in reverse transcription. A total of 3 1 plantlet samples were 
used for RNA extraction by the RNAprep Pure Plant Kit 
(polysaccharides & Poluphenalics-rich) (Tiangen, Beijing, China). 
The quality of all RNA samples was analyzed by the Synergy H 1 
Microplate Reader Multi-Mode (Bio-Tek, Winooski, VT, USA) 
with a 260/280 ratio from 1.9 to 2.1 and a 260/230 ratio from 2.0 
to 2.5. The integrity of the RNA samples was analyzed by agarose 
gel electrophoresis. 

The first-strand cDNA was synthesized with a 10 |iL reaction 
system according to the instructions of the TAKARA PrimeScrit 
RT Reagent Kit (Perfect for Real Time) (TaKaRa Biotechnology, 
Dalian, China). The quality and integrity of cDNA were 
determined in the same way as above. All the cDNA samples 
were diluted to 5 ng-uL 1 for the following qRT-PCR reaction 
and stored at — 20°C until use. 

qRT-PCR and Data Analyses 

An ABI 7500 Real Time PCR machine (Applied Biosystems, 
Foster City, CA, USA) and its default program (2 min at 50°C and 
10 min at 95°C followed by 40 cycles at 94°C for 15 s, and at 
60°C for 60 s.) were employed for qRT-PCR with a reaction 
mixture volume of 20 |jL in an optical 96-well plate. 10 (4.L of 
SYBR Green Master Mix (Roche), 10 pM of each primer, 10 ng 
of final cDNA and 6.4 |jL of RNase-free water were added to the 
reaction mixture. A control was also included in each plate with 
2.0 uL of RNase-free water as a template. Three technical 
replicates were contained in each plate. Specificity verification of 
the PCR amplification dissociation and the PCR efficiency curves 
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were determined for each candidate reference gene prior to the 
qRT-PCR evaluation of these genes in sugarcane. 

Statistical Analysis of Gene Expression and Comparison 
of Normalization Methods 

The PCR amplification efficiencies of all gene assays were 
estimated based on the slope of the standard curve in conjunction 
with the following formula: 

PCRefficiency (E) = lO^ 1 /*' 0 ^ - 1 (1) 

Standard, PCR efficiency and the correlation coefficient (R 2 ) 
curves of each gene were generated in Microsoft Excel 2003 using 
a range of seven dilutions made in ten, five or three-fold 
decrements on the YC05-179 control sample. GAPDH and TUB 
were ten-fold; 25S rRNA, ACT, 18S rRNA, eEF-la, CUL and eIF-4a 
were five-fold and CAC, TIPS-41, UBQ, APRT and PRR were 
three-fold. 

The overall expression levels of the candidate genes were 
transformed into threshold cycle (Ct) values by the ABI 7500 Fast 
Real Time PCR System. Ct values over 40 indicated undetectable 
product and were considered as missing values (NA) for 
subsequent calculations. After collecting and converting the 
threshold PCR cycle data, Ct average values for the evaluated 
genes were inputted into software according to the corresponding 
manuals of geNorm (trial version; Biogazelle, Zwijnaarde, 
Belgium) [2] and NormFinder (ver. 0.953) [1]. After that, the 
values of stability of candidate genes achieved from geNorm, 
NormFinder and deltaCt (only the deltaCt results from RefFinder 
were employed [16]) were used to calculate Pearson correlation 
values (r value) by SAS S21.0. The values reflect the level of 
correlation of between the results from geNorm, NormFinder and 
deltaCt. 

Results 

Screening of Candidate Reference Genes and Primer 
Design 

A total of 13 candidate reference genes were selected in 
sugarcane or in other plant species for evaluation on the basis of 
their stable expression across developmental stages and/or abiotic 
stresses. These included four previously assessed sugarcane 
candidate genes, 25S rRNA, GAPDH, ACT and TUB [31,32], as 
well as nine new candidate reference genes, APRT, PRR, 18S 
rRNA, eEF-la, eIF-4a, UBQ, CAC, TIPS-41 and CUL. Since there is 
limited available sugarcane genome sequence information, the 
publicly available gene sequences from O. sativa (18S rRNA, 
AK059783; eEF-la, AK061464; eIF-4a, AK073620; UBQ5, 
AK061988), Z mays {CUL, GRMZM2G166694_T04) and A. 
thaliana [CAC, TIPS-41) were used as the probes to search within a 
sugarcane expressed sequence tags (ESTs) database (www.ncbi. 
nlm.nih.gov/ nucest/?term = sugarcane). Two candidate reference 
genes in sugarcane which performed well in B.juncea [12], named 
CAC and TIPS-41, were identified by querying homologous 
sugarcane sequences with A. thaliana genes complete CDS [CAC, 
At5G46630; TIPS-41, At4G34270). The remaining five new 
candidates, 18S rRNA, eEF-la, eIF-4a, UBQ and CUL, performed 
well in 0. sativa [9] and Z ma J> s [H]- AU seven EST sequences, 
including those of 18S rRNA, eEF-la, eIF-4a, UBQ5, CAC, TIPS-41 
and CUL, were acquired from the publicly available database in 
NCBI using candidate ESTs with the highest homology to the 
target sequences. Both reference sequences from 0. sativa, Z mays, 



B. juncea, T. aestivum or A. thaliana and the corresponding target 
sequences from sugarcane were aligned together in DNAMAN to 
identify the complete sequence identity. This information was used 
to design primers using the Primer-BLAST tool from the NCBI 
(http://www.ncbi.nlm.nih.gov/tools/primer-blast/). The primer 
sequences for all 13 candidate reference genes are shown in 
Table 1. 

Verification of Primer Specificity, Efficiency and Gene 
Expression Profile 

Using the control sample of YC05-179 as the cDNA template, 
the specificity of primers used in qRT-PCR experiments was 
confirmed based on a melting curve analysis and agarose gel 
electrophoresis. The standard curve, PCR efficiency and the 
correlation coefficient (R 2 ) of each gene were generated in 
Microsoft Excel 2003 using a range of seven dilutions made in 
ten (GAPDH and TUB), five (25S rRNA, 18S rRNA, eEF-la, eIF-4a 
and CUL) or three-fold (CAC, TIPS-41, UBQ, APRT and PRR) 
decrements on the YC05-179 control sample. The qRT-PCR 
efficiency formula (Eq. 1) was used in the calculation. The 
regression coefficient R 2 for all the primers varied between 
0.9876-0.9999 over a serial dilution of cDNA. qRT-PCR 
efficiencies of primers ranged from 93.24% to 113.83% (Table 1). 

Over all samples, the Ct values of the 1 3 candidate reference 
genes varied over a wide range, and the mean Ct values of these 
genes across all the samples ranged from 14.27 to 28.21 (Table 1). 
Among these candidate reference genes, 25S rRNA was the most 
abundantly expressed gene in all of the samples (mean 
Ct±SD= 14.27±0.71) followed by 18S rRNA (mean 
Ct±SD= 15.38±0.74), whereas PRR was the least abundantly 
expressed gene (mean Ct = 28.21 ±2.04). ACT, GAPDH and eEF- 
la were close in Ct values (25.03, 24.64 and 24.24, respectively) 
but differed in SD values (1.06, 1.23 and 1.33, respectively). The 
mean Ct±SD of the remaining genes (TUB, TIPS-41, CAC, CUL, 
eIF-4a, UBQ, APRT and PRR) varied from 26.47 ±1.29 to 
28.21 ±2.04. The data also showed that the expression of ACT 
and eIF-4a were the least variable (Covariance (CV) of 4.23% and 
4.42%, respectively). The CV values of 18S rRNA, UBQ, 25S 
rRNA, GAPDH, CAC and ^PffTranged from 4.81% to 5.05%. The 
CV values of TUB, eEF-la and CUL ranged from 5.31 % to 5.52%. 
In contrast, the CV values of TIPS-41 and PRR were 7.67% and 
7.25%, respectively (Table 1). 

In order to make a wider expression analysis of candidate 
reference genes, the 57 diverse samples were divided into four 
experimental sets. The I s set consisted of leaf, leaf sheath, stem 
epidermis, stem pith and bud samples from "ROC20, 
"ROC"22, FN40, Liucheng03-182 and YC05-179 (Table 2). 
The 2 nd set was comprised of 12 samples from FN40, Liucheng03- 
182, "ROC'20 and YC05-179 treated with ABA, MeJA and SA. 
The 3 rd set contained 19 samples treated with H 2 O z , NaCl, PEG, 
CuCl 2 and CdCl 2 (abiotic stresses). In the 4 th set, both the 2 nd and 
3 rd sample sets were included. The variation of transcript levels 
indicated that the expression of the candidate genes was affected 
by tissue types and experimental conditions (Table 2). In the 1 st 
set, 18S rRNA performed as the least variable gene, whereas 25S 
rRNA, CUL and PRR had a larger range of expression level 
variation. 25S rRNA and eIF-4a had the lowest expression 
variation in the 2 nd set and, 25S rRNA and ACT had the lowest 
variation in the 4 th set. A CT and TUB performed less variability in 
the 3 rd set. Additionally, TIPS-41 &TUB and TIPS-41 &PRR 
displayed the most variable expression profiles when treated with 
hormones and abiotic stresses, respectively (Table 2). A significant 
variation in the expression of PRR and APRT was observed in all 
treatment samples (Table 2). Although none of the candidate 
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£ 3! reference genes displayed constant expression levels throughout 

g "<!f the various cultivars and among the different kinds of treatments, 

^ $ six genes among these 13 reference genes (25S rRNA, 18S rRNA, 

« 2 GAPDH, ACT, UBdand eIF-4a) varied at a relatively lower level 

?? J" according to CV values (Table 2). 
S"! 

| v Expression Stability Analysis and Ranking of All 13 

1 8 Candidate Reference Genes in Sugarcane 

g "g Based on the expression stability analysis of the 1 3 candidate 

rs reference genes by geNorm, the six top ranked genes in different 

[j 5 cultivar tissues were CUL= CAOAPRT>TIPS-41>GAPDH> 

g gi eEF-la. The top ranked genes in the cultured plantlets under 

o -fi hormone treatments were CUL= eEF- 1 a> GAPDH> GAG, and the 

Li = rank of the top six genes in the 3d set was CAC= CUL> APRT> 

? g GAPDH>eEF-la (Fig. 1). The top six genes that were expressed 

§ relatively constantly in the 4 th set were e£F-la= GAPDH>CUL> 

% Q CAC>APRT>eIF-4a (Fig. 1). It is interesting to note that the most 
stable traditional housekeeping gene was GAPDH compared to 25S 

£ rRNA, ^4CTand TUB, which is similar to previous reports [32,33]. 

Three reference genes, the eEF-la, GAG and CUL, were more 

E -5. stable than the other six new candidates (18S rRNA, eIF-4a, UBQ 
TIPS-41, APRT and PRR). Also, in response to various abiotic 
stresses and hormone treatments, UBQ_or TIPS-41 were identified 

2 5 as the least stable genes (Fig. 1). 

g ^ NormFinder analysis revealed that eEF-la and eIF-4tx were 

£j u always two of the top three most stable reference genes in the 

£ c treated samples, while the third place was occupied by ACT, 

I r5 GAPDH or TUB in the 2 nd , 3 rd and 4 th sets, respectively (Table 3). 

gj u Three reference genes, GAPDH, eEF-la and 18S rRNA, were most 

■is S stable in the different tissues from all the five tested cultivars 

-7. rj (Table 3). UBQ had the worst stability of the 13 candidate 

j" z reference genes in untreated samples and in samples under 

^ o hormone treatments (Table 3). Likewise, TIPS-41 ranked at the 

o „ bottom of these candidate genes in both the 3 rd and 4 th sets 
(Table 3). Though ACT, TUB and GAPDH expressed at less 

£ ?. variable levels in the 1 st , 2 nd and 4 th sets, the expression of all the 

8 t four traditional housekeeping genes fluctuated across all four 

.y g experimental sets and was less stable than eIF-4a or eEF-la in most 

>, -a of the samples (Table 3). 

"§ ™ The deltaCt method was used in RefFinder, which integrates 

■2 T geNorm, NormFinder, BestKeeper and deltaCt into a web-based 

"5 program. The ranking of the 13 genes based on the deltaCt was 

& _g consistent with that of NormFinder in different tissue samples 

.s except for the ranking of 25S rRNA (Table 3 and Table 4). 

y. o Comparing the average expression stability of the top six genes 

^ E valued by geNorm, NormFinder and deltaCt, the three candidate 

| °" genes GAPDH, eEF-la and eIF-4a performed better in all the 31 

o o treated samples, and eEF-la, APRT and GAPDH performed better 



a e 



S £ 



in all the 25 tissue samples (Table 4). Conversely, using the three 



g 2 algorithms, TIPS-41 and UBQwtrt the worst performing genes of 

~ S all treated samples and untreated samples, respectively (Table 4). 

•4= £ rs Pearson correlations were calculated among the geNorm, 

S -J °. NormFinder and deltaCt methods using the stability values in 

-o S § the 1 st and 4 th sets. Interestingly, the Pearson correlations among 

■f *t g all three stability tests were positive. A significant correlation was 

55 w observed between outcomes from NormFinder and deltaCt 

"= ™ 3 (r = 0.946), indicating that the ranking results of all 13 reference 

5^1 genes from the above two methods were nearly identical in the 1 st 

J s . ^. set (Table 5). The lower correlations (Table 5) between geNorm 

> I ft and NormFinder in the 4 th set (r = 0.438) and the 1 st set (r = 0.476) 

S3 5 were reflected by the results of the different rankings (Table 4). 
to 1/1 t: 
2? 0 

^ r\i "D 
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Figure 1. Gene expression stability of 13 candidate genes in sugarcane as predicted by geNorm. Average expression stability (M) 
following stepwise exclusion of the least stable gene across all the samples within an experimental set. The least stable gene is on the left, and the 
most stable on the right. The name elF-4a in the figure stands for elF-4x. ACT stands for "p-actin" and TUB stands for "p-tubulin". 
doi:1 0.1 371 /journal.pone.0097469.g001 



Optimal Number of Reference Genes for Expression 
Normalization across Different Experimental Sets 

The optimal number of reference genes in the normalizing 
experiment was determined by geNorm by calculating the pairwise 
variation (V n /V n+ i) between the normalization factors (NF) across 
all the samples of the different experimental sets. The pairwise 
variation (V = V n /V n+1 ), which was counted between NF n and 
NF n+ i, was used to find the best combination of genes for reliable 
normalization [2]. Vandesompele et al. suggested a threshold of 
V = 0.15, which indicates that adding one more gene has little 
influence on the calculation of normalization factor [2] . As shown 
in Figure 2, the rank order of gene stability established by means of 
stepwise exclusion of the least stable gene suggested that the 
combination of CAC, CUL, APRT and TIPS-41 could provide a 
dependable result while normalizing the qRT-PCR data of the 
target gene in the 1 st sample set (Fig. 2). However, in the 2 nd set 
the two reference genes [CUL and eEF-la) were enough to achieve 
a V 2 /3 = 0.134, which is close to the V = 0.15 threshold proposed 



by Vandeompele et al [2] and is thus an efficient and economical 
strategy to quantify sugarcane samples from hormone treatments. 
Similarly, the use of CAC and CUL was enough to achieve a V2/3 
= 0.144 in samples treated with compounds eliciting abiotic stress 
(Fig. 2). The V3/4 value in the 4 th set samples from four sugarcane 
cultivars was 0.101 which suggests that the combination of 
GAPDH, eEF-la and CUL is the best choice for quantification of 
gene expression in qRT-PCR (Fig. 2). 

Discussion 

Sugarcane, sorghum, maize and rice belong to the Andropo- 
goneae tribe and share many similarities in their genetic 
composition [36-40]. Rice, O. saliva, which diverged from a 
common Andropogoneae ancestor around 50-70 million years 
ago, performs as a model plant for grass species in several fields of 
modern molecular biology research [41]. Modern sugarcane 
cultivars are hybridized crosses between S. officinarum and S. 
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Table 3. Expression stability of 13 reference genes in Saccharum officinarum as calculated by Normfinder. 



Rank 


1 st set 




2 nd set 




3 rd set 




4 th set 




gene 


stability 


gene 


stability 


gene 


stability 


gene 


stability 


1 


GAPDH 


0.053 


elF-4a 


0.288 


elF-4x 


0.269 


elF-4ct 


0.218 


2 


eEF-la 


0.148 


ACT* 


0.385 


eEF-la 


0.462 


eEF-la 


0.419 


3 


1 8S rRNA 


0.156 


eEF-la 


0.466 


GAPDH 


0.506 


TUB** 


0.419 


4 


elF-4a 


0.167 


CUL 


0.557 


TUB** 


0.583 


UBQ 


0.445 


5 


APRT 


0.253 


PRR 


0.573 


UBQ 


0.620 


GAPDH 


0.505 


6 


25S rRNA 


0.255 


GAPDH 


0.583 


CUL 


0.666 


ACT* 


0.546 


7 


CAC 


0.256 


APRT 


0.615 


APRT 


0.671 


CUL 


0.571 


8 


TIPS-41 


0.258 


1 8S rRNA 


0.656 


ACT* 


0.680 


APRT 


0.580 


9 


CUL 


0.302 


CAC 


0.738 


CAC 


0.701 


PRR 


0.616 


10 


TUB** 


0.315 


25S rRNA 


0.816 


PRR 


0.729 


CAC 


0.668 


11 


ACT* 


0.335 


TIPS-41 


0.931 


25S rRNA 


0.748 


25S rRNA 


0.744 


12 


PRR 


0.417 


TUB** 


0.961 


18S rRNA 


0.918 


1 8S rRNA 


0.787 


13 


UBQ 


0.494 


UBQ 


0.990 


TIPS-41 


1.352 


TIPS-41 


1.092 



*P-actin, 
**fj-tubulin. 

doi:1 0.1 371 /journal.pone.0097469.t003 

spontaneum, resulting in a high degree of polyploidy and frequent 
aneuploidy [42]. Comparisons involving sugarcane as well as 0. 
sativa, S. bicolor and Z- mays indicated significant common 
conservation of gene content and a few rearrangements [40]. 
The results of the present study also showed that the stable 
expression of candidate genes, such as eIF-4cc, may follow the same 
pattern in S. officinarum, O. sativa, S. bicolor and £ mays. 

qRT-PCR is a quick, reliable and accurate tool to analyze gene 
expression. Accurate mRNA normalization requires one or two 
internal control genes which are stably expressed should be used. 



Since no single control is appropriate for all experimental 
treatments, it is generally suggested to select suitable internal 
controls prior to use for normalization. Although it is generally 
considered that housekeeping genes are expressed constantly, Lilly 
et al. and Nicot et al. revealed that the expression of such genes 
can undergo significant stability changes during biotic or abiotic 
stresses [18,43]. Thus, several statistical algorithms, such as 
geNorm [2], NormFinder [1], BestKeeper [17], the deltaCt 
method [10] and the RefFinder WEB-based software [11], were 
developed to assess expression stability. Among these algorithms, 
geNorm, NormFinder and the deltaCt method were used in the 



Table 4. Comparison of the expression stability of 13 reference genes in Saccharum officinarum as calculated by geNorm, 
Normfinder and deltaCt. 



Rank 


1 st set 












4 th set 












geNorm 




NormFinder 


deltaCt 




geNorm 




NormFinder 


deltaCt 




gene 


stability 


gene 


stability 


gene 


stability 


gene 


stability 


gene 


stability 


gene 


stability 


1 


CAC 


0.343 


GAPDH 


0.053 


GAPDH 


0.940 


GAPDH 


0.333 


elF-4a 


0.218 


CAC 


1.060 


2 


CUL 


0.343 


eEF-la 


0.148 


eEF-la 


1.040 


eEF-la 


0.333 


eEF-la 


0.419 


elF-4a 


1.070 


3 


APRT 


0.464 


18S rRNA 


0.156 


18S rRNA 


1.070 


CUL 


0.438 


TUB" 


0.419 


GAPDH 


1.080 


4 


TIPS-41 


0.657 


elF-4x 


0.167 


elF-4t>L 


1.100 


CAC 


0.453 


UBQ 


0.445 


CUL 


1.090 


5 


GAPDH 


0.726 


APRT 


0.253 


TIPS-41 


1.130 


APRT 


0.526 


GAPDH 


0.505 


eEF-la 


1.100 


6 


eEF-la 


0.777 


25S rRNA 


0.255 


APRT 


1.130 


elF-4a. 


0.598 


ACf 


0.546 


UBQ 


1.110 


7 


elF-4a. 


0.827 


CAC 


0.256 


CAC 


1.190 


PRR 


0.800 


CUL 


0.571 


APRT 


1.160 


8 


18S rRNA 


0.876 


TIPS-41 


0.258 


CUL 


1.220 


25S rRNA 


1.080 


APRT 


0.580 


25S rRNA 


1.220 


9 


TUB" 


0.955 


CUL 


0.302 


25S rRNA 


1.300 


ACf 


1.254 


PRR 


0.616 


18S rRNA 


1.390 


10 


25S rRNA 


1.023 


TUB" 


0.315 


TUB" 


1.300 


18S rRNA 


1.378 


CAC 


0.668 


ACf 


1.480 


11 


PRR 


1.102 


ACf 


0.335 


ACf 


1.430 


TUB" 


1.518 


25S rRNA 


0.744 


PRR 


1.930 


12 


ACT 


1.172 


PRR 


0.417 


PRR 


1.560 


UBQ 


1.602 


18S rRNA 


0.787 


TUB" 


2.190 


13 


UBQ 


1.274 


UBQ 


0.494 


UBQ 


1.810 


TIPS-41 


1.827 


TIPS-41 


1.092 


TIPS-41 


2.750 



*P-actin, 
**P-tubulin. 

doi:1 0.1 371 /journal.pone.0097469.t004 
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Table 5. Correlation coefficients based on the visualizing reference genes ranked by geNorm, Normfinder and deltaCt. 







1 st set 


4 th set 


geNorm VS NormFinder 


0.476** 


0.438** 


geNorm VS deltaCt 


0.679** 


0.674** 


NormFinder VS deltaCt 


0.946*** 


0.549** 



The p-value indicated by asterisks (***<0.0001; 
doi:1 0.1 371 /joumal.pone.0097469.t005 



*<0.05). 



present work. The deltaCt results were obtained from RefFinder, a 
web-based tool that integrates the most commonly used statistical 
algorithms. The normalization results were mosdy affected by 
variation of the quantity and quality of introducing RNA and 
cDNA or even reaction-to-reaction variations [1]. 

When employing a variable reference gene, Gutierrez et al. 
revealed that nearly 100-fold variations could be found when 
quantifying target gene expression, thus leading to a huge potential 
scope for misinterpretation of the expression pattern of target 
genes [7]. Guo et al. observed much different expression values 
during a study of sugarcane dirigent protein gene expressed in 
sugarcane stems when using the different internal control genes of 
25S rRNA and GAPDH [44]. Therefore, to achieve a reliable result 
with qRT-PCR, a systematically validated reference gene should 
be used and taken as an essential component of qRT-PCR 
analysis. In the present study, the three statistical algorithms, 
geNorm, Normfinder and deltaCt, were used and the Ct values 
were inputted into RefFinder. When correlation analysis based on 
the ranking order of the 13 candidate genes was applied for 
comparison of the three statistical algorithms, the results showed 
that the correlation coefficient between NormFinder and deltaCt 
was more positive (significantly in the cultivar set, r = 0.946) than 
the correlations between geNorm and Normfinder or between 
geNorm and deltaCt. These results were in accordance with the 
results of Jacob et al. [45]. The ranking order of the evaluated 



candidate genes by NormFinder and deltaCt in our study was 
generally consistent. 

According to the CV values, 25S rRNA, 18S rRNA, GAPDH, 
ACT, CAC, eIF-4a, UBQand APRT showed low variation in gene 
expression in all four experimental sets (Table 2). However, only 
GAPDH, eEF-la and eIF-4a performed well in all three statistical 
algorithms. These three genes had relatively stable expression in 
different sugarcane genotypes and tissues. Our results also showed 
that the three genes CUL, eEF-la and elF-4a were suited for gene 
expression normalization in hormone-treated experiments, while 
GAPDH and eEF-la were ideal when analyzing samples compa- 
rable to the 3 rd sample set. However, when these two sets (the 2 nd 
and 3 rd sets) were integrated into the 4 th set and analyzed together, 
eEF-la, eIF-4oc and GAPDH were the most stable reference genes. 

In previous studies, GAPDH was found to possess the most stable 
mRNA expression in sugarcane [32] and 25S rRNA was the most 
stable gene in sugarcane infected with Ustilago scitaminea [33]. 
Similarly, this study identified GAPDH as one of the most stable 
reference genes for all types of hormone and abiotic treatments. 
25S rRNA and 18S rRNA were two of the more abundantly 
expressed target genes under investigation, which is contrary to the 
reference gene selection principle that emphasizes moderate 
expression levels. It is difficult to detect the variation of 25S rRNA 
and 18S rRNA due to their rich content, so these two genes may 
not be ideal for expression normalization. It should be pointed out 
that some genes with stable but low expression levels, such as genes 



0.25 






V2/3 




V3/4 




V4/5 




V5/6 




V6/7 


1 1 


V7/8 




V8/9 




V9/10 




V10/11 




V11/12 




V12/13 



Figure 2. The optimal number of reference genes required for effective normalization in each of four experimental sets in 
sugarcane. The pairwise variation (V n /V n+1 ) was analyzed between normalization factors NF n and NF n+1 by geNorm program to determined the 
optimal number of reference genes for accurate normalization in samples from different sugarcane cultivar samples (1 st set), sugarcane hormone- 
treated (2 nd set), abiotic-treated (3 rd set) and treatments (hormone-& abiotic-treated, 4 th ). ACT stand for "P-actin" and TUB stand for "p-tubulin". 
doi:1 0.1 371 /journal.pone.0097469.g002 
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coding for certain transcription factors may be used as internal 
references in qRT-PCR experiments. A suitable reference should 
have a minimum difference in Ct value since this would lessen 
influence on the quantification and hence be more accurate [8,54]. 
In the present study, four historically used housekeeping genes, 
GAPDH, 25S rRNA, ACT and TUB, were also included in the 
evaluation [32,33]. 25S rRNA, ACT and TUB had poor 
performance across all four experimental sets of sugarcane 
samples, and similar poor performance was seen when using 
some of the novel candidates, such as UBQ TIPS-41, PRE and 18S 
rRNA. 

In this study, eEF-la, which has been proven to be a suitable 
reference gene for expression normalization in 0. sativa and C. 
sativus [5,12], also ranked at the top when evaluated by geNorm, 
NormFinder and deltaCt across our four experimental sets. Being 
an important specific protein factor involved in the process of 
protein translation, eIF-4oc had the same performance as eEF-la in 
most experimental treatment conditions in our study. Our study 
showed that eIF-4cc had a relatively high expression under most 
experimental conditions, which is in agreement with the study of 
Zhu et al. in C. papaya [46]. eIF-4cc also performed well in Musa 
paradisiaca, Lycoris longituba, Hevea brasiliensis and Cqffea spp. [47-50]. 
The low copy number of APRT 'm the sugarcane genome, which 
was reported by Casu et al. [30], provides an advantage in 
analyzing gene copy number. However, in the present study the 
expression of APRTwsls found to be easily affected by abiotic stress 
conditions (3 rd set) and thus its application is limited. With more 
copies in the sugarcane genome than APRT [30], PRR displayed 
variable performance under different stresses and in different 
sugarcane tissues. 

In other plants such as B.juncea [15], S. lycopersicum [22], P. edulis 
[29], C. papaya [46], G. max [51] andZ. culinaris [52], the expression 
profiles of housekeeping genes, such as GAPDH, eIF-4a and eEF-la 
were not expressed as consistently as non-traditional housekeeping 
genes such as CAC, UBQ9, TIPS-41, MTB, ELFlb/dOs, TBP1, 
TBP2, TIF, RPL2 and PP2Acs. However, in the present study 
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